Hybrid Monte Carlo Simulation of Graphene on the Hexagonal Lattice 
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We present a method for direct hybrid Monte Carlo simulation of graphene on the hexagonal 
lattice. We compare the results of the simulation with exact results for a unit hexagonal cell 
system, where the Hamiltonian can be solved analytically. 
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INTRODUCTION 

Graphene, a single layer of carbon atoms forming a 
hexagonal lattice, has remarkable properties [T]. In the 
tight-binding approximation the quadratic Hamiltonian 
gives origin to a dispersion formula which for low mo- 
menta is analogous to the dispersion formula for relativis- 
tic fermions in two dimensions [2]. This has prompted 
some researchers to adapt to graphene lattice gauge the- 
ory techniques which have been profitably used for the 
study of Quantum Chromodynamics and other parti- 
cle systems [3]. In the work of [4] one approximates 
first the tight-binding Hamiltonian of graphene with a 
Dirac Hamiltonian, incorporates the Coulomb interac- 
tion through the introduction of a suitable electromag- 
netic field, and finally discretizes the resulting continuum 
quantum field theory on a hypercubic space-time lattice. 
The hybrid Monte Carlo method [5], widely used in lat- 
tice gauge theory to simulate fermions interacting with 
quantum gauge fields, can then be used to investigate 
the effects of the Coulomb interaction in the graphene 
system. 

The approach outlined above has led to very interest- 
ing and valuable results [4], yet one would think that, 
since the starting point is a system already defined on a 
lattice, it should be possible to apply the hybrid Monte 
Carlo technique directly to the graphene lattice. The 
clear advantage of this approach is the direct connection 
to the experimentally determined physical lattice con- 
stants of the tight-binding model, which represents an 
accurate description of the experimental system. In this 
letter we illustrate how this can be done. 

Graphene is a system of interacting electrons located 
at the vertices of a hexagonal lattice. It is convenient to 
think of the graphene lattice as consisting of two trian- 
gular sublattices, which we denote by A and B, which 
together with the centers of the hexagons (sublattice C) 
form a finer, underlying triangular lattice (Fig. [I]). We 
introduce fermionic annihilation and creation operators 
a XyS ,a x s for the electrons on the two sublattices, where 
a; is a site index and s — ±1 is the spin index. The lat- 



FIG. 1: The hexagonal graphene lattice consists of the two tri- 
angular sublattices A (solid) and B (empty), which together 
with the centers of the hexagons (dots) form a finer, underly- 
ing triangular lattice. 



tice must be made finite in order to perform numerical 
simulations. While there is a broad range of boundary 
conditions of physical interest, here we consider periodic 
systems formed by identifying opposing sides of a hexag- 
onal lattice of length L, illustrated in Fig. [I] for L = 4. 

The tight-binding Hamiltonian H consists of two 
terms: the quadratic kinetic term, 



+ h.c.) 



(1) 



where the sum runs over all pairs (x, y) of nearest neigh- 
bor sites (coupling the A and B sublattices) and the two 
values of the spin, and the Coulomb interaction Hamil- 
tonian, 



where 



He = ^e 2 V Xiy q x q y , 



1x — a x l a x,l + a x -l a x,-l 
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(2) 



(3) 



is a local charge operator and V is the interaction poten- 
tial. We have explicitly introduced the charge coupling 
constant e. 
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Several comments are in order. First note that in 
the kinetic term we have neglected the smaller next-to- 
nearest neighbor hopping within each sublattice, which 
would introduce a small (probably manageable) complex 
phase in the path integral. The charge operator Eq.[3]has 
a —1 to account for the background charge of the carbon 
ion: it ensures that the system is neutral at half filling, 
and it will play an important role for our functional in- 
tegral formulation. V could be the actual 3d Coulomb 
potential, but could be any other interaction potential. 
The only thing crucial for us is that the matrix V x , v be 
positive definite. Finally, we note that the Hamiltonian 
of Eqs. T][2 commutes with the isospin generators 



a x , s >/2. (4) 



In order to explore the properties of the system one 
would like to calculate expectation values, 



(Ol(tl)0 2 (*2 



Z- 1 TrT[0 1 (t 1 )0 2 (t 2 )e- 



(5) 



where j3 can be interpreted as an extent in Euclidean 
time, T[. ..] stands for time ordering of the opera- 
tors inside the square bracket with respect to the Eu- 
clidean evolution implemented by exp(— f3H), and Z = 
Tr exp(— fiH) is the partition function. 



PATH INTEGRAL FORM 

Our goal is to provide an equivalent path integral for- 
mulation of Eq. [5] conducive to calculation by numeri- 
cal simulation, following a rather standard procedure to 
convert from the Hamiltonian into a Lagrangian. Wc 
will first express the expectation values and the parti- 
tion function in terms of an integral over anticommut- 
ing fermionic fields, i.e. elements of a Grassmann alge- 
bra. (The literature on the path integral formulation of 
quantum expectation values is very rich. In our work we 
followed the very clear and useful formulation given in 
the first chapters of [§].) This gives origin to an inte- 
grand with an exponential containing a quadratic form 
in the fermionic fields, from H2 and the normal ordering 
of He, as well as a quartic expression, from He- The 
quartic expression can be reduced to a quadratic form 
by a Hubbard-Stratonovich transformation [7], through 
the introduction of a suitable auxiliary bosonic field (in 
our case a real field) , and now the Gaussian integral over 
the elements of the fermionic variables can be explicitly 
performed, leaving an integral over the bosonic field only. 
The problem, however, is to obtain an integral that can 
be interpreted as an integration over a well defined prob- 
abilistic measure, which can thus be approximated by 
stochastic simulation techniques. We will show here how 
the symmetries of the system make this possible. 

We start by rewriting the expression for the charge as 



We now introduce hole creation and annihilation opera- 
tors for the electrons with spin — 1: 



b\ = a x _!, b x 
so that the charge becomes 



l x,—l 



q x = a x a x 



(7) 



(8) 



Note that we dropped the spin indices since from now on 
a, and b, W will always refer to spin 1 and —1, respec- 
tively. Finally we change the sign of the b, w operators 
on one of the sublattices. The crucial constraint is that 
all redefinitions of the operators respect the anticommu- 
tator algebra. From the fact that H2 only couples sites 
on the two different sublattices, it follows that H2 now 
takes the form 



H 2 = ^ + b t b v + h - c 0- 

We introduce fermionic coherent states 

e -E a; (>/' : c<4+ I fe b x)|o) j 



(9) 



IVS»7) 
(ip*,r]*\ 



(10) 



where ip x ,ip*,r) x ,r) x are anticommuting fermionic vari- 
ables (elements of a Grassmann algebra). 

The path integral formulation is obtained by factoring 

e -p H = e -H5 e -HS £ -H6 (jy t terms) (11) 

with 5 = f3/N t , and then inserting repeatedly among the 
factors the resolution of the identity expressed in terms 
of an integral over the fermionic variables. The trace in 
Eq. [5] must also be expressed in terms of a similar inte- 
gral. (See e.g. [6] for details.) This leads to integrals over 
fermionic fields t/j x j,ip* t ,rj x t , 77* t (the index t = 0, N t — 1 
appears because of the multiple resolutions of the iden- 
tity and can be thought of as an index labeling Euclidean 
time), which contain in the integrand expressions of the 
type 



The last ingredient is the identity 

x,t> Vx,t) 



(12) 



(13) 



which is true of any normal ordered function F of the 
operators a\,,b\,,a x ,b x . 

The Hamiltonian is in fact already in normal order 
except for the local term e 2 V xx q x q x , which can be written 
as the sum of two normal-ordered pieces, 



,i a x,i — a Xt -\a x 



(6) 



*V xx q x q x = e 2 V xx : q x q x : + e 2 V xx (a x a x + b\,b x ), (14) 
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By reassigning the quadratic term in Eq. 14 to H2, the 
exponent — H 5 in Eq. [12] is normal ordered but the ex- 
ponential exp(— H 8) is not. However exp(— H S) differs 
from its normal ordered form : exp(— H 5) : by terms 
0{8 2 ). So, in the limit of N t ~ > 00 one may replace the 
operator expression exp(— H 8) with an exponential in- 



volving the fermionic fields, as follows from Eq. 13 This 
leads to the following expression for the partition func- 
tion 



,N t -l 



(15) 



Z = lim / TT dipmdtpmdri^drjm 

N t ->oo J 

m— 

x e -E m .„(^M m ,„^+, ? ,*„M m ,,^„) e -E ;c , ! ,, t e 2 Qx, t K c , !( Q B , tI 5 

where Q xt = ip x tip x ,t ~ Vx tVx.t and we have used m (or 
n) as a shorthand for the indices x,t. M is a matrix 
whose components may be deduced from 



+ e 2 V xx ip* it ip X:t - k Y2 (ipl, t 4>v,t + ^*y,ti>x,t) S 



where ip Xi N t m ust be identified with —ipx,o- 

We now perform a Hubbard-Stratonovich transforma- 
tion, introducing c-number real variables cf> Xt t to recast 
the exponential with the quartic term in the form 



e > Qv,tV x , v Q v , t 6 



n« 



(17) 



where we have absorbed a constant measure factor in the 
definition of the integral over (^re- 
inserting the r.h.s. of Eq. 17 into Eq. 15 we get 

J x,t 



lx,tU<Px,t 



(18) 



It is convenient to introduce a matrix $ which is diagonal, 
with diagonal entries 

$*,t = <fc,t& (19) 
With this, Eq. [18] can be written in very compact form 



Z = I dtfidip* dipdrj* dr/ 



x e 



-4>V 1 (pS/i-xp* (M+te$)ip-r)* (M— je*)?j 



(20) 



where we have used matrix notation for all the sums and 
have dropped the limit notation. 



The Gaussian integration over the anticommuting vari- 
ables can now be done to obtain 



Z 



d(f>e 



-(f>V~ <j>6/4 



det(M-ze$)det(M + ie$). (21) 



Because of the identity, 

det(M-ze$)det(M + «e$) = det[(M + ze$) f (Af + ie$)] 

the measure is positive definite. The down spins are 
treated as antiparticles (holes) moving backward in time 
relative to the up spins, exactly canceling the phase for 
each separately. Correlators for the fermion operators 
are now be obtained by integrating the appropriate ma- 
trix elements of (M + ie$) _1 or (M — ze<I>) -1 with the 
measure given by Eq. |21[ 

Equation [21] is the main result of our work. It estab- 
lishes the partition function and expectation values as 
integrals over real variables with a positive definite mea- 
sure. This is a crucial step for the application of stochas- 
tic approximation methods. There remains the problem 
of sampling the field 4>x,t with a measure which contains 
the determinant of a large matrix. But, following what is 
done in lattice gauge theory, this challenge can be over- 
come through the application of the hybrid Monte Carlo 
(HMC) technique [5]. In a broad outline, in HMC one 



first replaces the determinants in Eq. 21 with a Gaussian 
integral over complex pseudofermionic variables C, x ,t- 

det[(M + «e$) f (M + ie$)] 

= f dCd(e~ C{M+2e ® }i '^(A'+^r'C _ ( 22 ) 



(In this equation and in the following Eq. 23 we absorb 
an irrelevant, constant measure factor in the definition 
of the integrals.) One then introduces real "momentum 



variables" Tr x ,t conjugate to 4> Xj t and inserts in Eq. 22 
unity written as a Gaussian integral over tt. One finally 
arrives at 



Z = / d<f>dTrdCd( 



x e 



-<j)V- 1 4>8/A-C(M+ie<5>) ] _1 (M+ie*) _1 C-Tr 2 /2 



(23) 



The idea of HMC is to consider the simultaneous distri- 
bution of the variables 4>,ir,C and C* determined by the 
measure in Eq. [23j The phase space of these variables 
is explored by first extracting the ir, £ and £* according 
to their Gaussian measure, and then evolving the <fi and 
7r variables with fixed £, according to the evolution 
determined by the Hamiltonian 



+ C(M + ie$y- L (M + ie<i>y L C- 

(24) 

Because of Liouville's theorem, the combined motion 
through phase space will produce an ensemble of vari- 



ables distributed according to the measure in Eq. 23 and, 
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in particular, of fields <p distributed according to the mea- 



sure of Eq. 21 



Of course, the discussion above assumes that the 
Hamiltonian evolution of and ir is exact, which will 
not be the case with a numerical evolution. The HMC 
algorithm addresses this shortcoming by: 1) approximat- 
ing the evolution with a symplectic integrator which is 
reversible and preserves phase space, 2) performing a 
Metropolis accept-reject step at the end of the evolution, 
based on the variation of the value of the Hamiltonian. 



NUMERICAL TESTS 

We tested our method on the two-site system obtained 
by taking L = 1, which can be solved exactly. We label 
the sites x — 0, 1. With k = 1/3, the Hamiltonian H = 
H-2 + He is now 



H 2 = -(a\a + ajai + b\b Q + 6q&i) 
+ (i(%ao + blb + a\a\ + b\bi) 

H c = 2e 2 (4a - blbo)(a\ ai - b\b{) ■+ 



(25) 



2e2 

TO 



where we have taken Vo i = Vj. o = 1/3 and a local in- 
teraction term Vo,o = Vi,i = 1/ro- The radius ro sets 
the physical scale in lattice units for localization of the 
net charge at the carbon atom. It must be restricted 
to tq < 1 for stability of the vacuum. Also the nor- 



mal ordering prescription for e V xx q x q x in Eq. 14 adds 
a new contribution to Hi in the form of an /3 "chemi- 
cal potential" /j,a x s a| s a X:S '. It is well known [5] that an 
23 chemical potential for any value of /i does not intro- 
duce a phase in the measure. To maintain the full SU (2) 
"flavor" symmetry of the tight-binding graphene Hamil- 
tonian, we must set ^ — e 2 /r a . For the two-site system, 
the isospin generators of Eq. [4] become 



(-l) a a!>t and h = ^ + 6 t 6 J/ 2 - l 



allowing us to unambiguously classify the 16 states as de- 
generate I = 0, 1/2, 1 isomultiplets: 5 singlets, 4 doublets 
and one triplet. 

We compared HMC results for expectation values of 
several products of fermionic operators with the corre- 
sponding exact values, finding satisfactory agreement. 
For example, the correlation function 



C a (t) = ((00 - oi)(t) (oj - 4)(0)>/2 



(26) 



is illustrated in Fig. [2j which shows HMC results con- 
verging to the exact correlators for both the free theory 
with e = and an interacting case with e = 0.5. 

A stringent test is to demonstrate the convergence to 
exact SU(2) symmetry in the "time" continuum limit. 
To this end, consider a second correlation function, 

C b (t) = ((bl + b\)(t)(b + b 1 )(0))/2, (27) 



o.ui 
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FIG. 2: HMC data with N t = 64 and 128 converging to exact 
results for the correlation function of Eq. 26 with /1 = e 2 /ro, 



ro = 1/2 and fixed /3 = N t S — 6.4. The upper set of results 
are for e = 0, the lower set for e = 0.5. 




0.01 



0.001 



FIG. 3: HMC data with N t = 64 and N t = 128 for the 



correlation functions of Eqs. 26 (""["", red squares and blue 



triangles) and |27| ("j", pink squares and green triangles) con- 
verging to the same exact continuum result (black solid line) . 
Simulation parameters are the same as in Fig. [2] 



related to C a (t) by an SU(2) rotation. In Fig. [3] we com- 
pare HMC and exact results for both C a (t) and Cb(t), 
finding that both correlators converge to the same con- 
tinuum result. 

We extract the energies of the isodoublet states at 
nonzero i5 = P/N t by fitting the correlator data in 
Fig [3] to single exponentials, C(t) s» eT Et for fit range 
0.4 < t < 4. The results in Fig. [4] clearly show linear be- 
havior E ~ Eq + ci<5, converging to the exact continuum 
energy is E Q = e 2 + y/4 + e 4 - 1 « 1.266. The con- 
tinuum limit is consistent with restoration of the SU(2) 
symmetry of the Hamiltonian: a simultaneous linear fit 
to both sets of energies gives lim^o E = 1.262 ± 0.004, 
with ci = 1.25 ±0.07 and -0.98 ±0.04 for correlators C a 
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FIG. 4: Linear extrapolation (including error band) of the 
HMC energies for isodoublet correlators C a it) (blue dots) and 
Cb (red triangles) as a function of the "time" lattice spacing 
8 = 6.4/iV t for N t = 32, 64, 128 and 256. The dotted hori- 
zontal line marks the exact continuum result. 



modated by reweighting without a prohibitive cost. The 
crucial observation is the cancellation between the phase 
of the up spin and down spin determinant, when the lat- 
ter are treated as holes moving backward in time. We are 
currently pursuing simulations of larger systems, and be- 
ginning to explore the many possible investigations and 
generalizations (e.g., distortions of the lattice, phonons, 
inclusion of magnetic fields). 
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and C{,, respectively. 

CONCLUSION 

While the results we reported are for a small test 
system, they demonstrate that HMC simulations of 
graphene directly on the hexagonal graphene lattice are 
possible and have the potential to produce valuable re- 
sults. The dominant nearest neighbor hopping term has 
no sign problem, and we anticipate that a small next- 
to-nearest neighbor coupling k'/k ~ 0.05 can be accom- 
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